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We study average flow properties in porous media using a two-dimensional network simulator. It 
models the dynamics of two-phase immiscible bulk flow where film flow can be neglected. The 
boundary conditions are biperiodic which provide a means of studying steady state flow where com- 
plex bubble dynamics dominate the flow picture. We find fractional flow curves and corresponding 
pressure curves for different capillary numbers. In particular, we study the case of the two phases 
having equal viscosity. In this case we find that the derivative of the fractional flow with respect 
to saturation is related to the global pressure drop. This result can also be expressed in terms of 
relative permeabilities or mobilities, resulting in an equation tying together the mobilities of the two 
phases. 

PACS number(s); 47.55.Mh 



I. INTRODUCTION 

Two-phase displacements in porous media have dur- 
ing the past decades been studied by experimental work 



statistical models 
The work in this 



numerical simulations 
-UJ, and differential equations 
field has been motivated by applications in oil recovery, 
and also hydrology. The complex nature of transport pro- 
cesses in porous media makes it worthwhile to approach 
the problem in many ways. 

In this paper we present simulation results based on a 
network simulator of a porous medium jl^ . The model 
is on the pore level scale in the sense that each pore is 
represented in the network. The generic building blocks 
are tubes which constitute all inter-pore connections and 
which hold the pore-space volume. As will become clear 
in section each tube allow for a maximum of two 
fluid-fluid interfaces. This captures essential properties 
of porous media throats. These aspects set the level of 
coarse graining of the system. There has been done nu- 
merous work on lower length scales and mesoscale meth- 
ods, such as lattice gas methods An overview over 



some different scale models can be found in 15 1. On 
the other hand, compared to large scale reservoir simu- 
lations our work is detailed and it may provide av- 
erage properties which can be used as input parameters 
in work on larger scales. Currently, results are mainly 
generic, and more detailed work is required to produce 
results which are valid for specific porous media. In par- 
ticular, simulations in 3D are required. 

In general, the model can be used for 3D porous media 
and irregular node positions. This follows from the fact 
that each inter-pore throat of the real porous medium 
is replaced by a tube in our model. However, in this 
paper we restrict ourselves to the study of flow in 2D 
using a regular network with random properties. Pre- 
vious use of the model in 2D is simulation of drainage 
The previous simulations have successfully re- 
produced the temporal evolution of the three regimes in 
drainage: viscous flngering, capillary fingering, and sta- 



ble displacement. The major important innovation in the 
present work is the change of boundary conditions. Be- 
fore there was an inlet and an outlet between which an 
invasion process was simulated. Now biperiodic bound- 
ary conditions are used, which makes the system closed, 
and the system is run for a longer time. This gives av- 
erage statistical fiow properties in a new way, and it is 
the biperiodic boundary conditions which are the major 
difference compared to other somewhat similar network 
models @,|, see Sec. 

The nature of the simulations is so that the systems 
approach a steady state, in which the two phases flow in a 
bubble mixture. The notions of drainage and imbibition 
are not adequate to describe the flow. On pore level 
there is a continuous process of breaking up bubbles and 
merging together bubbles. Experimental work showing 
similar complex bubble dynamics have been done in two- 
dimensional etched glass networks [p^l-pl]. 

The basic results which are presented are fractional 
flow curves of the nonwetting phase as a function of non- 
wetting saturation. Also we present the corresponding 
global pressure drops which are applied. This is done for 
six different capillary numbers. The fractional flow and 
pressure drop can be transformed into the terminology 
of relative permeabilities and mobilities. 

Our simulations show that for the case where the two 
phases have equal viscosity, there is a relation between 
fractional flow and pressure. This relation can be ex- 
pressed in the form of a first order differential equation; 
the derivative of the fractional flow with respect to satu- 
ration is connected to the global pressure drop, also as a 
function of saturation. 

The paper is organized as follows. Section || provides a 
brief description of the model. Section |l| concerns sim- 
ulation results, starting out with d escribing the nature 
of the simulations. Thereafter III A de als w ith fractional 
flow and pressure curves in general. III B discuss how 



these curves depend on the capillary number, and III C 
establish the relationship between the derivative of the 
fractional flow and the the global pressure drop. Finally, 
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in [V we make concluding remarks. 



II. MODEL 

The basis for the present work is a network model 
which represents the two-dimensional porous medium. 
The network model is based upon the one presented by 
Aker et al. [ pTjp^ . The geometric representation of the 
porous medium and the possible spatial fluid distribution 
is essentially unchanged. Some changes have been made 
though. They concern the boundary conditions and the 
detailed rules for motion of interfaces. All the details 
concerning these poi nts have been presented in a previ- 
ous piece of work flSt] . Therefore we will restrict ourselves 
here to an essential resume of the physically important 
aspects and omit the computational details. 

The porous medium is represented by a square lattice 
of cylindrical tubes, tilted 45° with respect to the overall 
flow direction. The volume of throats and pores is con- 
tained in these tubes. The points where four tubes meet 
are referred to as nodes. Randomness is incorporated in 
the system by allowing the position of each node to be 
randomly chosen within the interval plus minus thirty 
percent of the lattice constant away from its respective 
lattice point. From these positions the distance dij be- 
tween connected nodes i and j is calculated for all i and 
j. Further, the average radius of each tube is chosen by 
random in the interval (0.1c?, 0.1d-|-0.3dy), where d is the 
average distance between nodes, i.e. the lattice constant. 
Clearly, other topologies in 2D (or in 3D) can be chosen. 
However, as a first approach to the study of the relation- 
ship between other variables it is necessary to work with 
one fixed topology. We will address the topology ques- 
tion again in the discussion. There current restriction to 
2D is also very convenient because 3D systems are CPU 
time demanding. 

We consider two fluids within this system of tubes. 
They are separated by a set of interfaces, menisci, in the 
tubes. We do not allow for film flow. Motion of the fluid 
during a simulation is represented by the motion of the 
menisci. In each tube we allow zero, one or two menisci. 
If at any instant the evolution of the system generates a 
third meniscus in one tube, then the three menisci are 
collapsed into one. The position of this new meniscus is 
the one which preserves the volumes of the two fluids. 
This upper limit of two menisci in each tube sets the 
resolution of the fluid distribution. 

With respect to permeability the tubes are treated as 
cylindrical, but with respect to capillary pressure they 
are hour-glass shaped. This means that over each of the 
menisci in the system there is a capillary pressure which 
varies with the menisci's position in the tubes. The for- 
mula for the capillary pressure is 



which is a modified form of the Young-Laplace]^ law 
p5| , ^ . Here, r is the radius of the tube, and 7 is the 
inter-facial tension between the fiuids. Further, x is the 
position of the meniscus in the tube, running from zero 
to one. 

The volumetric fiux through one single tube is given 
by the Washburn equation for capillary flow [E3|; 
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where pc is given by Eq. (|l|). Considering the tubes as 
cylindrical with radius r, the fraction irr^ / d is the cross- 
sectional area divided by the length of the tube. The 
permeability of the tube is fc = which is known from 
Hagen-PoiseuUe flow. When two fluids are present in a 
tube, an effective viscosity /ioS is used. The viscosities 
of the two fluids are weighted according to their volume 
fraction within the tube at the beginning of the time 
step to give /ieff- Ap is the pressure difference between 
the ends of the tube. 

In a network of tubes, the net flux passing through the 
nodes must be zero. There is assigned a pressure to each 
node. Using the Washburn equation (|^) for the flow in a 
single tube, the flux conservation gives a large system of 
linear equations in the pressure variables. The equations 
must be solved using the desired boundary conditions 
which we will describe shortly. 

The traditional use of network simulators has been the 
study of invasion processes. Usually the flrst row of nodes 
is an inlet row where all the nodes have the same flxed 
pressure. Similarly the last row works as an outlet row, 
typically having all nodes fixed at zero pressure. For a 
given time step the pressures in these two rows are fixed 
constraints. The pressures in all the other nodes are free 
variables that are solved for. 

Once the pressure field is known. Equation (^ gives 
the flux in each tube. In turn this information is used 
to forward integrate the system one time step, using the 
explicit Euler scheme. In practice, forward integration 
means motion of the menisci, and as long as the menisci 
move within a single tube this is straightforward. How- 
ever, there will be menisci reaching the ends of the tubes. 
This is dealt with as follows for every node separately. All 
the menisci which have proceeded past the edges of the 
tubes neighboring a node, and so have entered the node, 
represent incoming volume of one of the fluids. In order 
to have volume conservation the same amount of volume 
must leave the node. This is done by adding new menisci 
in the neighboring tubes in which the flow is away from 
the node. The positions of these new menisci are such 
that volume is preserved. This is the essential part of 
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*Young-Laplace law is sometimes referred to as Laplace's 
equation psl]. 
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these rules of motion. However, the rules which are ac- 
tually used are a little more complicated, mainly due to 
computational technicalities. They have been presented 
in detail before jl3|, and they are not necessary for the 
understanding of results in this paper. 

Many authors have studied invasion under a constant 
applied pressure. In the work by Aker et al. constant 
flux invasion has been considered. By solving the flow 
field for two different globally applied pressures giving 
two different fluxes, one can calculate the pressure that 
would give the desired flux. We do not go into details 
here since this is described in detail before ||l^. We will 
just remark for now that all simulations presented in this 
paper are done with constant flux rate. 

Network models based on Washburn's equation was 
introduced by the Payatakes group, see e.g. and ref- 
erences therein. Their work related to steady-state two- 
phase flow is of particular interest [^,^ . A detailed com- 
parison of the two lines of modeling is not our aim, but 
we give some general remarks. Our model is a little more 
coarse in two ways. First, the resolution of pore space and 
of the fluid in the pore space is lower in our work since it 
is a tube model, and since only two menisci are allowed in 
each tube. Second, film flow is not included, neither ex- 
plicit as network components or via the use of coefficients 
for, e.g., coalescence of oil blobs. Also, we have chosen to 
work with fixed values for a larger number of parameters. 
In other words, we present generic rather than specific use 
of the model in this piece of work. Finally, we empha- 
size that our model use biperiodic boundary conditions 
to reach steady-state flow. A short description follows. 



A. Boundary conditions 

Invasion simulations go on until the invading fluid 
reaches the outlet. If they were to go on further, they 
change character since one fluid is percolating the sys- 
tem. We wish to address the question of what happens 
in a system far from inlets and outlets, i.e. given a very 
large system, we take out a small piece somewhere in the 
middle and study its properties. 

This is done by adjoining the inlet row and outlet row 
so that the fluid that flows out of the last row enters 
the first row. In practice this works in such a way that 
the simulation can go on for ever, regardless whether one 
fluid percolates the system or not. In a sense having 
biperiodic boundary conditions make the system infinite. 
However, the system is closed, and there is a fixed volume 
of each of the fluids in the system. Thus each simulation 
takes place at constant saturation equal to the one of 
the initial configuration. We return to this point when 
presenting results from simulations in section |l|. 

These boundary conditions can be visualized by con- 
sidering the system to be the surface of a torus. The flow 
is driven around the torus by a numerical trick. This 
method has been used in connection with random resis- 



tor networks ||2^. Briefly, one can say that there is an 
invisible line dividing the system. Looking over this line, 
there is a pressure fall or jump, which is the driving force 
of the flow. All the tedious details are given in the pre- 
vious piece of work none of which are necessary to 
know in order to understand the contents of the present 
work. 



III. SIMULATIONS 

The essence of the model porous network is that it 
is situated on the surface of a torus. In other words 
the applied boundary conditions are biperiodic, which 
makes the system closed. By means of the pressure fall 
technique, the overall flow is around the torus, while the 
total volumes of the two fluids remain constant within 
the system. 

Equally important is the fact that the simulations, 
which are presented here, are done by keeping the to- 
tal flux around the torus constant. As time evolves the 
fluid distribution changes and local capillary pressures 
change. In order to keep the total flux constant, the 
globally applied pressure fall must change in every time 
step. A sample of global pressure versus time is given in 
Fig. g(a). 
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FIG. 1. The figures show the typical time evolution of (a) 
global pressure, and (b) flux of nonwetting fluid, in a system 
of size 20 X 40 nodes. The nonwetting saturation is approx- 
imately 35%, the capillary number is Ca ~ 3.2 x 10~^, and 
the total flux through the system is Qtot ~ 14.0 x 10~^dm^/s. 
On the average the fluid travels almost 17 times around the 
system during this simulation. 

We wish to draw attention to two aspects of the nature 
of this curve. The first is that it appears noisy. On a short 
time scale local redistributions of the fluids and changes 
in the capillary pressures lead to a change in the overall 
pressure which must be applied in order to produce a 
given flux. Also one should remember that the time axis 
is very compressed. The fluctuations in pressure are not 
achieved in one single time step, but rather in the order 
of one hundred time steps. 

The second aspect is that on a longer time scale the 
systems have a transient part and a steady part. The 
transient part depends on the initial conditions, while in 
general we have found that the steady part is independent 
of the initial conditions. After the systems reach the 
steady state, the properties of this state are typically the 
time averages over the variables: global pressure, the flux 
of each of the fluids, the number of interfaces, velocity 
distributions, and others. 

Fig. 0(b) shows the nonwetting flux Q function 
of time. The data in |l|(a) and |l|(b) are from the same 
simulation. We see how the nonwetting flux also has a 
transient and a steady part. More convenient is the non- 
wetting fractional flow, which is defined as the flux of the 
nonwetting fluid through the system divided by the total 
flux: Fnw = Qnw/Qtot- To some extent the fractional 
flow property has been studied before. Previous results 
will be summarized shortly. 



A. Fractional flow and pressure curves 



sion, 7 = 30.0mN/m. The cross-sectional area is approx- 
imately S — 0.145cm^ for the system size which is used. 
We have further chosen the value of the viscosity to be 
/i — O.lPa s. Actually the total flux is the only one of 
these parameters which we vary. This is no limitation as 
we will discuss shortly. 

For the capillary number, Ca — 3.2 x 10^"^ the nonwet- 
ting fractional flow as a function of nonwetting saturation 
is given in Fig. 0(a). The corresponding global pressure 
is shown in Fig. ^(b). Here it should be noted again that 
the values of the fractional flow and global pressure are 
the time averages in the steady state. 

The very important question of how these results de- 
pend on the system size was examined before Q . When 
all other parameters were held constant, except for sys- 
tem size and saturation, the fractional flow as a function 
of saturation showed to be independent of the system 
size. The examined system sizes were 20 x 40 nodes, 
20 X 80, and 40 x 80. 

The interpretation of the size independence is that the 
system is large enough to be in the asymptotic limit. Of 
course, this might happen to be untrue for some possible 
parameter sets of the system, but for the sets used in 
the present work it is true. All simulations which are 
presented here are therefore done with system size 20 x 40 
nodes. 




There are many possible parameters of the model sys- 
tem. We will generally use the saturation as the indepen- 
dent variable in our plots. The nonwetting saturation is 
deflned as the volume of nonwetting fluid divided by the 
total volume, Snw — Kiw/Vtot- The other possible pa- 
rameters are the viscosities of the fluids, the total flux, 
the inter-facial tension, and system size. 

In this work we have restricted ourselves to the case 
where the two fluids have equal viscosity. The viscosity 
and some other parameters are combined in the capillary 
number 
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where ^ is the viscosity, Qtot is the total flux, 7 is the 
inter-facial tension, and E is the cross-sectional area of 
the system. Physically this number is the ratio between 
typical viscous forces and capillary forces within the sys- 
tem. We have used a fixed value for the inter-facial ten- 
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FIG. 2. The plots show average values of (a) fractional 
nonwetting flow and (b) global pressure, where the aver- 
ages are taken in the steady part. The capillary number is 
Ca = 3.2 X 10"^ We see clearly how these values depend on 
the nonwetting saturation. The simulations are run on five 
different geometries generated by different random seeds. 

The system size dependency of the global pressure is 
slightly more complicated. The pressure is applied in 
such a way that there is a pressure gradient along one 
of the system axes. Any change in size sideways, that is 
perpendicular to the gradient, does not alter the pressure 
curves. However, increasing the length along the pressure 
gradient with a factor results in a global pressure which 
is increased by the same factor. This is just another 
way of saying that it is the gradient in global pressure, 
which must be applied to give a desired total flux, that 
is independent of system size. In this work we consider 
only one system size so we prefer to use AP for the global 
pressure drop. 

The points in the plots come from five different geome- 
tries. They are drawn from the same distributions with 
different random seeds. The specific permeability of the 
porous medium, denoted by k, is related to the applied 
pressure and the total flux for a single phase, Darcy's 
law, 

Q^kAP^^ (4) 
E ji L ' 

where L is the length of the system along the pressure 
gradient, and the subscript on APg denotes single phase 
pressure. This is a geometrical property and will there- 
fore vary slightly for the five geometries. In fact this 
means that the pressure curves in Fig. ^(b) should be 
five slightly different curves. The variations in the spe- 
cific permeability are of the order of a couple of percents. 
The five curves are scaled according to this difference. 
The average position of the curves is kept constant un- 
der the scaling. We note that this implies that the sin- 
gle phase points, 0% and 100% saturation, are shown as 
single points. The rest of the curves are not perfectly 
overlapping due to the statistical variations. 

Subsequently we will take this scaling one step further. 
For illustrative purposes the physical dimension of pres- 
sure was used in Fig. ^(b). However, it will be more con- 
venient to work with a normalized pressure. This is done 
by giving the pressure in units of the single phase pres- 
sure, i.e., letting the axis be AP/APg- This is straighfor- 
ward, but note that the single phase pressure is a function 
of the capillary number. It follows from Eqs. (H) and (^) 
that APs cx Ca. 

The fact that the capillary number, as it is defined in 
Eq.(||), serves as the relevant combined parameter was 
also established in the previous piece of work This 
means that one can make changes to all the variables in 
Eq. (^) , but as long as the capillary number is preserved 
the shape of the fractional fiow curve will be the same. 

The characteristics of the fractional flow curve are as 



follows. Below a certain nonwetting saturation, approxi- 
mately 5% for Ca — 3.2 X 10^"^, we observe that the non- 
wetting fractional flow is essentially zero, meaning that 
only the wetting phase flows. Conversely above a certain 
nonwetting saturation, 85% for Ca — 3.2 x 10"'^, only 
the nonwetting phase flows. Between these values both 
phases flow. 

Imagine there had been no capillary forces between the 
two phases. Then both phases would have flown equally 
easy through the porous network. The nonwetting frac- 
tional flow would then be exactly equal to the nonwet- 
ting saturation. This is illustrated by the diagonal line 
in Fig. ^(a). The effect of having interfaces with capil- 
lary pressures can be thought of as the change from the 
ideal diagonal to the obtained curve. For low nonwetting 
saturations the fractional nonwetting flow is lower than 
the saturation, i.e. under the diagonal. For high nonwet- 
ting saturation the nonwetting fractional flow is higher 
than the saturation, i.e. above the diagonal. Roughly we 
can say that the phase of which there is more volume, 
gains and flows more than one naively might expect. At 
some point the curve cross the diagonal, and this is the 
point where neither phase gains compared to its volume 
fraction. This point is not at 50% saturation, and it is 
clear that there is an asymmetry between the wetting 
and nonwetting phases. 

The pressure curve in Fig. |2|(b) has two points plotted 
with 'x', at 0% and 100% saturation. These are single 
phase pressures, and since the two fluids have the same 
viscosity the pressures are equal. For all other satura- 
tions, two phases are present within the system. There 
are interfaces between the phases with capillary pres- 
sures. Motion of two fluids with interfaces through tubes 
requires more global pressure than for the case of a single 
fluid. One might imagine that in some cases some or all 
of the present interfaces do not move, which is typically 
the situation when only one fluid flows. Still the global 
pressure must be higher than for a single fluid because 
the nonmoving interfaces block out tubes, reduce possi- 
ble pathways, and hence reduce the specific permeability 
of the medium. As we can see from the figure, the global 
pressure is higher in the entire two-phase region than in 
the single-phase points. 

The general nature of the pressure curve is that with 
increasing saturation the pressure increases monotoni- 
cally to a maximum value. Thereafter it decreases mono- 
tonically. The curve in Fig. ||(b) is generally smooth 
except at nonwetting saturation of approximately 85% 
where the pressure drops rather rapidly. We observe for 
now that this point is coinciding with the point were the 
nonwetting fractional flow becomes unity. Also the char- 
acter of the curve changes at this point. To the left it is 
curved while to the right it is nearly a straight line. 
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B. The dependency on Ca 

Now we have learned that the capillary number serves 
as the relevant parameter for the fractional flow versus 
saturation curves. For two fluids having equal viscosity 
we have performed simulations for six different capillary 
numbers, which is shown in Figs. and ^. The system 
size is 20 x 40 in all simulations. This size is sufficiently 
large to be in the asymptotic limit, but still so small 
that simulations can be done within reasonable amounts 
of time. For each capillary number we have chosen 19 
different saturations from approximately 5% up to 95% 
in steps of 5%. Further we have used five random seeds 
which provide five different geometries. 

The shape of the fractional flow curve depends strongly 
on the capillary number. For high capillary numbers the 
curve is approaching a straight line. Intermediate capil- 
lary numbers, as for the sample which was discussed in 
the previous subsection, the curve has a clear S-shape. 
For lower capillary numbers the curve seems to approach 
a step function. 

The curves can be divided into three regions. The first 
being for low nonwetting saturation where the nonwet- 
ting fractional flow is zero. For the two highest capillary 
numbers this region is almost vanishing. As the capillary 
number is lowered, this region grows in size. The third 
region is for high nonwetting saturations where the non- 
wetting fractional flow becomes unity. The same general 
behavior is valid here as for the first region. It becomes 
wider with decreasing capillary number. The difference 
is that for a fixed capillary number the region with unity 
nonwetting fractional flow is larger than the region with 
zero nonwetting fractional flow. This means that there 
is not perfect symmetry between the two phases. 

The second region is the central part where both phases 
are mobilized. The width of this region decreases with 
decreasing capillary number. In the limit of inflnite cap- 
illary number, negligible capillary forces, the curve will 
approach the diagonal and thus span all of the saturation 
range. In the opposite limit of small capillary number, 
it is reasonable to expect that the curve will approach 
a step function, although we have not explicitly checked 
this for lower capillary numbers than presented here. Be- 
tween these limits there is a range of capillary numbers 
where the central parts of the curves have interesting 
structure. 
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FIG. 3. The figures show the nonwetting fractional flow as 
function of nonwetting saturation for six difi'erent capillary 
numbers. For high capillary numbers the curve is close to 
the diagonal, while for low capillary numbers the curve is 
S-shaped. The range of saturations for which both phases are 
mobilized decreases with decreasing capillary number. 

As it was pointed out in the previous subsection for 
Ca = 3.2 X 10^"^, the curves lie above the diagonal for 
large nonwetting saturations and below for small nonwet- 
ting saturations. The interpretation is that the system 
favors transport of the phase of which there is more vol- 
ume. This is not exactly true, since the crossover from 
favoring one phase to the other is not at 50% nonwetting 
saturation, but at a somewhat smaller value. Again the 
system is not perfectly symmetric in the two phases. Ac- 
tually the crossover point is a function of the capillary 
number. For very high Ca the crossover point is close 
to zero. It approaches 50% when the capillary number 
decreases. This situation corresponds to the problem of 
bond percolation in two dimensions, for which the per- 
colation threshold is known to be Sc = 1/2 |2^. In per- 
colation theory Sc is a critical bond probability, and it 
corresponds to a critical saturation in our problem. For 
small systems, history may evolve in such a way that the 
systems have a saturation other than 50% even though 
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only one phase flows. However, on the average, or in an 
infinite system which is self-averaging, the limiting value 
is 50%. We have listed the values of the crossover point 
in Table | 

Corresponding to the fractional flow curves in Fig. ^ 
are the global pressure curves in Fig. ^ Also here the 
curves depend strongly on the capillary number. It was 
pointed out in the previous subsection for d, = 3.2 x 10"'^ 
that the curve increases monotonically to a maximum 
value and thereafter decreases monotonically. This is 
true for the entire range of capillary number which we 
have examined. The position of the maximum value in- 
creases with decreasing capillary number. The values are 
listed in Table |. 

The pressure curves can be divided into three regions 
just like the fractional flow curves. At least for the three 
lowest capillary numbers, Fig. ^(d)-(f), their boundaries 
are clear and distinguishable. The first region is the one 
were the pressure increases linearly with nonwetting sat- 
uration from the value at zero nonwetting saturation. For 
the three highest capillary numbers. Fig. ^(a)-(c), this 
region is almost vanishing. 

The third region is where the curves decrease linearly 
with increasing nonwetting saturation towards the value 
of pressure at unity nonwetting saturation. This region 
can be seen in the entire range of capillary numbers. By 
inspection we see that the region boundaries are the same 
in the pressure curves as for the fractional flow curves. 
Thus the same comments regarding the width of the re- 
gions are valid for the pressure curves. 

The second region is the central part of the curve. It 
is curved and has nontrivial structure. For the lower 
capillary numbers the pressure increases abruptly at the 
boundaries of the region. In this region both phases are 
mobilized. The immediate conclusion is that more pres- 
sure is needed when both phases are mobilized, as in- 
terfaces then are mobilized. The relation between the 
fractional flow and the global pressure in this region will 
be the subject of the next subsection. 

The flrst and the third region share several properties. 
At their outer boundaries, 0% and 100% respectively, the 
values of the global pressure are equal, and equal to the 
single phase pressure. Also at the inner boundaries of 
these regions the global pressures have the same value. 
That is to say, as far as both inner boundaries are distin- 
guishable, we can make this observation. Physically the 
inner boundaries are just at the saturation values where 
there is an onset of mobilization of the other phase. The 
fact that the pressures are the same here is the same as 
saying that in this respect there is symmetry between the 
two phases. 
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FIG. 4. The figures show the global pressure, '-I-', as a 
function of nonwetting saturation for six different capillary 
numbers. The global pressures are normalized with respect 
to the single phase pressure APs as defined in (^. Recall 
that for each capillary number the total flux is fixed to a 
constant value. The curves correspond to the fractional flow 
curves in Fig. ^, and the derivate of the fractional flow curves 
are included here in a scaled version, 'o'. This scaling and 
the r elationship between the curves are treated in subsection 
[II Cl 



Within these regions one of the phases is immobilized. 
Still there is a linear dependency of the global pressure 
on the saturation within the regions. The reason for this 
change in pressure cannot come from the increased pres- 
sure necessary to move an increased number of interfaces 
since the interfaces generally are pinned to fixed positions 
when only one phase flows. When there is an increased 
volume of the phase which is immobilized, then there is 
less available volume for the other phase to move in. This 
is a geometrical constraint. The specific permeability of 
the single phase which flows will decrease within these 
two regions as the saturation of other phase increases. 
Thus more and more pressure is needed to maintain a 
constant total flux. 

This qualitative explanation does not account for the 
fact that these parts of the curves are straight lines. If 
the immobilized phase were pinned to positions in the 
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network which were more or less random, then to a first 
approximation a hnear increase in the saturation of the 
immobile phase would lead to a linear decrease in the 
effective system size for the mobile phase. In turn that 
would lead to a linear increase in global pressure. 

The next question is why the third region is wider than 
the first region. According to the first approximation 
reasoning above, they should be equal. Going a little 
beyond this approximation, we know that a bubble of 
nonwetting fluid is more likely to get pinned around a 
node. Its most stable position is when it is bounded by 
menisci in all neighboring tubes, and when they are in 
the half parts of the tubes which are closer to the node. 
Whatever direction the bubble is pushed it will have to 
pass the threshold in the respective tube. Bubbles of 
wetting fluid have opposite preference. Their most sta- 
ble position is within a tube. With two bounding menisci 
placed on each side of the center of the tube, the wetting 
bubble is stable to fluctuations in both directions. Of 
course, bubbles come in a range of sizes and this asym- 
metry between the two phases is more pronounced for 
smaller bubbles. On the average we can say that the 
nonwetting phase when it is immobilized will block more 
nodes than tubes. The wetting phase when immobilized 
will occupy comparatively more tubes than nodes. It is 
more effective to block nodes than tubes, and this ex- 
plains why the nonwetting phase is more effective than 
the wetting phase in reducing the specific permeability 
of the other fluid. The result is that the pressure curve 
is steeper in the flrst region than in the third, and thus 
the flrst region is narrower than the third. 



C. Relationship between Fnw and AP 

This subsection is devoted to the study of the central 
range of saturation where both phases are mobilized. In 
this region we have defined the crossover point as the 
point where the fractional flow is equal to the saturation. 
On both sides of this point the phase whose saturation 
is higher gains in the sense that the fractional flow is 
higher than the saturation. In this region the fractional 
flow curves are roughly S-shaped. In particular, for the 
capillary numbers in Fig. ^(c)-(d), the curvature changes 
approximately at the crossover point. Whether this is 
coincidental will be discussed later. We have determined 
the derivative of all the six curves. From these data we 
have estimated the turning points of zero curvature which 
are listed in Table |. 

The differentiation has been done in the straightfor- 
ward way, except at the boundaries of the central re- 
gion. Outside the region the derivative is zero, and inside 
the region the curves are smooth. However, right at the 
boundary the fractional flow curves have, for the lower 
capillary numbers, a break where the derivative almost 
diverges. We have ignored these breaking points. 

We find that the derivative has the same shape as the 



globalpressure in the central region. This is shown in 
Fig. H where the derivative is marked by 'o'. By the 
same shape we mean that we for each capillary number 
can find two dimensionless constants A and B for which 



B(Ca) = 



AP 
AR 



(5) 



The quality of the overlap of the two sets of data is very 
good for the four highest capillary numbers. Fig. ^(a)- 
(d). For the fifth. Fig. |(e) the quality is fair. For 
the lowest capillary number the points start to become 
spread out so much that a comparison of shapes is dif- 
ficult. We have included the derivative here to get an 
estimate for the scaling coefficients A and B. The values 
of the scaling coefficients are listed in Table ||. 

In order to understand the meaning of Eq. (^, it is 
useful to rewrite the expression in the terminology of mo- 
bilities and relative permeabilities [^^. The nonwet- 
ting relative permeability fcr,nw is defined by 



X " 



where the fraction 



(6) 



(7) 



is the nonwetting mobility. Here the constant k is the 
specific permeability as it was defined in Eq. (^. Sub- 
stituting for fc, keeping in mind that both phases have the 
same viscosity /x, and expressing Q^w in terms of Fnwi we 
obtain 

APs 

^r,nw('5'nw) — -^nw('5'nw) X . , ^ r. (8) 

i'-'nw j 

Likewise the result for wetting relative permeability is 



^r,w(5'nw) — -^w(*S'nw) ^ 



APiSn 



(9) 



where the single phase pressure drop APg is the same be- 
cause we consider two fluids with equal viscosities. The 
relative permeabilities for the capillary number Ca = 
3.2 X 10"'^ are shown in Fig. ||. These curves are com- 
parable to experimentally obtained relative permeability 
curves which are frequently used in the petroleum Indus- 
try M. 
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FIG. 5. The nonwetting and wetting relative permeabil- 
ity is shown as a function of nonwetting saturation. They 
are defined in Eqs. (^) and (^. The capillary number is 
Ca = 3.2 X 10~^, and these curves correspond directly to the 
fractional flow and pressure curves which are found in Fig. ^ 

Analogously to nonwetting mobility in Eq. we 
define wetting mobility as 



nw ) 

and the total mobility as 



Af(S'nw) — Afnw(5'nw) + M^{Sn 



(10) 



(11) 



In the simulations the total flux has been held constant. 
It can be expressed by using Eq. and its wetting 

counterpart, as 



Qtot = Qnw + Qw = S(A/nw + M^)k X 

= X Af(5„w)AF(5nw), 



AP 
17 



(12) 



which can be solved for AP; in combination with (|j) this 
is right hand side of Eq. (||). Likewise the fractional flow 
can be expressed in terms of the mobilities, 



F — 
nw — 



M 



(13) 



Thus the scaling relation in Eq. (ra) can be rewritten as 



A X 



Afnw(*S'n 

M{S nw 



B = 



1 



^iM{Sn 



where A and B are independent of the saturation, 
ferentiating once more, we find 



(14) 



Dif- 



A X 



1 dM(5n„) 

X rr: -■ (15) 



dS„ 



The new insight given by this equation is as follows. Gen- 
erally the nonwetting and wetting mobility are considered 
as two independent functions of saturation. The result 
in Eq. (|l^) shows that the two mobilities are related 
through an equation. 

The dependence which we have found is so far re- 
stricted to the case where both phases have equal vis- 
cosity. Whether the result will be extendable to the case 
of two different viscosities in some form is an interesting 
question, but has been outside the scope of the present 
work. The validity is also restricted to capillary numbers 
in the range 3.2 x 10"'' < Ca < 3.2 x lO'^. For higher 
capillary numbers there is little reason to expect any new 
interesting behavior since the asymptotic behavior of the 
fractional flow curve is the straight diagonal. Lower cap- 
illary numbers than this range are more interesting, also 
from a practical point of view since they may occur in 



reservoir conditions. Challenges in this region of parame- 
ter space are increased history dependence of the results, 
and considerably increased CPU time. 

The validity of Eq. (^ comes from visual inspection 
of the data collapse in Fig. ^. Visually it seems that for 
each capillary number, the two curves have their maxi- 
mum value at the same saturation. We wish to discuss 
whether these maximum points are coincidental with the 
crossover points on the fractional flow curves. Estimates 
for all three points for each capillary number are listed 
in Tabic | The error estimates come from the data anal- 
ysis, and for the turning point precision is lost in the 
differentiation process. The crossover points are sensi- 
tive to systematic errors, in particular for high capillary 
numbers. The reason is that the fractional flow curve 
becomes increasingly parallel to diagonal for increasing 
capillary numbers. A small vertical shift of the frac- 
tional flow curve will give a large shift in the crossing 
point. From the values in the table we observe that the 
three points coincide for the four lowest capillary num- 
bers within the error bars given. That is, the turning 
point for the two lowest capillary numbers are very un- 
certain and not listed respectively. For the two highest 
capillary numbers the crossover point is not equal to the 
maximum pressure within the error bars, but again one 
may speculate in the possibility of systematic errors in 
these two points. 

In conclusion, we cannot say whether the crossover 
point is related with the maximum point. However, it is 
an interesting question because it is suggestive to look for 
effective theoretical descriptions starting at a point hav- 
ing this supposedly nice behavior. Experimental work on 
this problem will provide means of checking both Eq. (|^) 
and the possible relation between the crossover and the 
maximum pressure point. 

TABLE I. The table lists the crossover points of the frac- 
tional fiow curves (crossover), the maximum points of the pres- 
sure curves (maximum AP), and the turning points(turn) of 
the fractional flow curves for six different capillary numbers. 
The points are all saturation values, which are dimensionless 
numbers in the range from zero to one. The three points 
seem to coincide for the four lowest capillary numbers. The 
crossover point differs from the other two for the two highest 
capillary numbers. 



Ca 




crossover 


maximum AP 


turn 


3.2 X 10" 


-2 


0.13(2) 


0.19(4) 


0.22(4) 


1.0 X 10" 


-2 


0.18(2) 


0.22(5) 


0.23(6) 


3.2 X 10" 


-3 


0.25(2) 


0.28(5) 


0.25(7) 


1.0 X 10" 


-3 


0.33(2) 


0.31(4) 


0.34(4) 


3.2 X 10" 


-4 


0.40(2) 


0.38(5) 


0.31(8) 


1.0 X 10" 


-4 


0.48(3) 


0.45(6) 
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It is interesting to look at how the crossover point and 
scahng factors vary with the capillary number. We ar- 
gued that the crossover point should approach the perco- 
lation threshold Sc in the limit of small capillary numbers. 
This limit turns out to be not zero but a lower critical 
capillary number Ca,crit- We denote the crossover point, 
the crossover saturation, s. The crossover point as a func- 
tion of the logarithm of the capillary number is shown 
in Fig. 1^. That is, the capillary number is normalized 
with respect to the critical capillary number which is de- 
termined from fitting the following relation to the data 
points, 



a In 



(16) 



Here Sc is known so that the value of a = 0.066 is the 
slope of a straight line fit, and the value of Ca,crit = 
7.3 X 10^^ is the one which makes the line meet the s- 
axis at Sc — 0.5. 
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FIG. 6. The crossover points s from Table | are shown as 
a function of the logarithm of the capillary number. The 
straight line fit is made on the basis of the four lowest capillary 
numbers, which are closest to the critical capillary number. 

The physical idea behind having a critical capillary 
number is that at some point the viscous forces of the 
flowing fluids and the pressure gradient will become so 
small that they cannot mobilize any more interfaces. In 
that situation one of the phases will have a continuous 
pathway to flow in. The other phase is pinned to its 
current locations. Therefore any further decrease in the 
capillary number below the critical value, will not add 
anything to the picture; one phase flows. 

The scaling coefficients A and B from Eq. (^) are 
listed in Table ||. The dimension of these coefficients 
are the same as for pressure. The immediate tentative 
interpretation of the meaning of the two are as follows. 
The constant B is the threshold pressure which is applied 
right at the borders of the central region of the fractional 
flow curves, the pressure where there is onset of mobility 
of both phases. The factor A is the one that must be 



multiplied to the derivative of the fractional flow curve 
in order to obtain the pressure jump from the onset level 
B to the actual pressure at a given saturation in the 
central region. 
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FIG. 7. The scaling factor A from Eq. (||). Data are taken 
from Table 0. 
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FIG. 8. The scaling factor B from Eq. (^). Data are taken 
from Table The value of Be is the one which gives the best 
power law fit; Be ~ 0.65. 

It is clear that the pressure is an increasing function of 
the capillary number. Therefore it is natural that A and 
B are functions of the capillary number. The functional 
forms are presented in Figs. and |[ We can say that 
the explanation of these two power laws is related to the 
width of the central region. Also the extent to which 
interfaces are put in motion plays a role. Both these 
quantities depend on the capillary number. We have so 
far not been able to construct sound explanations for 
these power laws, and present them as observations only. 

TABLE II. The scahng factors A and B from Eq. (||) are 
shown for the six capillary numbers. They are used in the 
scaling of the derivative of the fractional flow curves which is 
shown in Fig. 0. 



Ca 



A(Ca) 



B(Ca 



3.2 X 10- 



0.22 



0.87 
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1.0 X 10"^ 0.35 1.21 

3.2 X 10"^ 0.83 1.84 

1.0 X 10^^ 1.47 3.83 

3.2 X 10"'' 3.22 6.90 

1.0 X 10"'' 5.89 17.6 



IV. DISCUSSION 

The problem of two-phase flow in porous media is com- 
plex. We attack the problem with a network simulator 
on pore level. The model is coarse grained on the level 
of the internal structure of the pores. The current use 
of the model is generic, but specific use is possible with 
few modifications. Useful average properties of porous 
media can be obtained which may be used as parameters 
on larger reservoir scale simulations. 

In this paper we present results for 2D square networks 
of tubes. This topology is chosen for convenience. It is 
possible to extend the work to 3D and choose both regu- 
lar and irregular connectivity. The choice of a particular 
topology is one aspect of making the model specific to 
a given porous medium. Further, average tube radius 
and the the width of the tube radius distribution are im- 
portant parameters that should be adjusted for the same 
purpose. Here, we have used a distribution which corre- 
sponds to experimental sizes in Hele-Shaw cells with glass 
beads which have been used experimentally . Further 
work and comparisons with these experiments may pro- 
vide insight into how the model can be callibrated in 
order to become quantitavely precise. 

The most important characteristic of the simulations 
in this paper, which sets them apart from other simula- 
tions, is the biperiodic boundary conditions. This makes 
the system closed with a fixed saturation of each of the 
two phases. For six different capillary numbers we run 
the systems until they reach a steady state where the 
flow is characterized by complex bubble dynamics. The 
notion of imbibition and drainage are not adequate to 
describe this situation, which would also be the situation 
deep inside reservoirs. We flnd average flow properties 
as a function of the saturation. These properties are the 
fractional flow and the global pressure drop, from which 
in turn one can calculate relative permeabilities and mo- 
bilities. 

In particular for the case of two phases having equal 
viscosity, we flnd that the derivative of the fractional flow 
is related to the global pressure drop, see Eq. (^. This 
relation ties together these two properties which gener- 
ally are believed to be independent. We have shown how 
this statement can be transformed into the language of 
relative permeabilities and mobilities, which as a conse- 
quence are not independent but must obey Eq. (p^. 

An important note here is that this result so far has 
only been established by numerical work. It is very in- 
teresting to get an experimental veriflcation of this re- 
sult. Hopefully, after an experimental verification, this 
equation can be of assistance in the measurement of two- 
phase flow properties. In the experimental situation it 
is difficult to measure all variables precisely. The satu- 
ration can e.g. be reconstructed from the pressure and 
fractional flow relationship and Eq. (||). 
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